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For open quantum systems coupled to a thermal bath at inverse temperature /3, it is well known 
that under the Born-, Markov-, and secular approximations the system density matrix will approach 
the thermal Gibbs state with the bath inverse temperature /?. We generalize this to systems where 
there exists a conserved quantity (e.g., the total particle number), where for a bath characterized 
by inverse temperature j3 and chemical potential n we find equilibration of both temperature and 
chemical potential. For couplings to multiple baths held at different temperatures and different 
chemical potentials, we identify a class of systems that equilibrates according to a single hypothetical 
average but in general non-thermal bath, which may be exploited to generate desired non-thermal 
states. Under special circumstances the stationary state may be again be described by a unique 
Boltzmann factor. These results are illustrated by several examples. 
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Thcrmalization is a classical phenomenon: Coupling 
two materials at different temperature will lead to equi- 
libration at some intermediate temperature - depending 
on the heat capacities of the constituents. Especially 
when one piece is significantly larger than the other, the 
temperature of the larger piece will hardly change, such 
that it may be understood as a heat bath. In contrast, 
the temperature of the smaller piece will simply approach 
the bath temperature in that limit. 

The dynamics of open quantum systems that are cou- 
pled to a thermal bath is however more difficile [1, 2]. A 
powerful tool to describe the evolution of such systems 
in various limits is the quantum master equation [3, 4]: 
A first order differential equation - typically with con- 
stant coefficients - describing the evolution of the sys- 
tem part of the density matrix. As the derivation of an 
exact master equation is impossible in most cases, one 
has to rely on perturbative schemes. In such schemes, it 
is often already a challenge to preserve the fundamental 
properties of the density matrix such as its trace, its self- 
adjointness, and its positive semidefinitencss. Starting 
from microscopic models, especially the last property is 
often hard to fulfill, as for master equations with con- 
stant coefficients, preservation of positivity requires the 
dissipator to be of Lindblad [5] form. Such Lindblad 
form dissipators are gcncrically derived in the singular 
coupling limit [6], the weak-coupling limit - also termed 
Born-Markov-secular [7] (BMS) approximation - and in 
coarse-graining schemes [8]. Within the BMS approxi- 
mation, thermalization of the system and equilibration 
of the systems temperature with that of the bath have 
been proven [9]. However, some baths are not only de- 
scribed by a temperature, but may also equilibrate under 
further side constraints - typically modeled by a chem- 
ical potential. When we consider couplings to multiple 
baths held at different temperatures [7, 10] and/or dif- 
ferent chemical potentials [11, 12], we have the generic 
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situation for transport [13] from one reservoir through 
the system to another reservoir, which may be used to 
generate interesting non-equilibrium stationary states in 
the system. 

This paper is organized as follows: Having introduced 
the terminology in Sec. I we show how conserved quanti- 
ties lead to additional properties of the dampening coef- 
ficients in Sec. II. The case of a single bath is discussed 
in Sec. Ill A, followed by a discussion of multiple baths 
in Sec. IIIB. We derive general statements on the result- 
ing non-equilibrium stationary state for master equations 
that are tridiagonal in the system energy eigenbasis in 
Sec. IIIC. The conditions under which such a stationary 
state may still appear thermal are discussed in Sec. Ill D. 
Finally, the results are demonstrated with a number of 
examples in Sec. IV. 

I. PRELIMINARIES 

We will consider a large closed quantum system with 
the total Hamiltonian 

H = H S + H B + H SB , (1) 

where H§ and H B act only on the system and bath parts, 
respectively, and -ffse mediates a coupling. The latter 
may generally be decomposed as [3] 

M 

H S b = Y. Aa ® Ba ( 2 ) 

a=l 

with M hermitian system (A a = A^ a ) and bath (B a = 
BJA coupling operators. By convention, the system cou- 
pling operators may be chosen traceless and orthonormal 
Trl^cA^} = 5 a p. For example, for an ./V-dimensional 
system Hilbcrt space one may use the M = N 2 — 1 gen- 
erators of the symmetry group SU(7V) for the system 
coupling operators. 

Under the Born, Markov, and secular (BMS) approx- 
imations [3], and assuming that the bath is kept in an 
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equilibrium state p B with the properties Ttb {-BqPb} = 
as well as [Hb,Pb] = 0, one derives a master equation of 
Lindblad form for the system density matrix p$. In the 
system energy eigenbasis H$ \a) = E a \a) it assumes the 
form [14] 
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where d a b = defines the unitary action of decoher- 
ence (also denoted Lamb-shift [3] or exchange field [15]) 
and the dampening coefficients "/ a b,cd describe the non- 
unitary (dissipative) terms due to the interaction with 
the reservoir. The net effect of the secular approxima- 
tion mentioned above is that these coefficients may van- 
ish when some transition frequencies are not matched 
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which is formally expressed by the Kronecker-5 symbols 
(but see [14, 16] for a Lindblad-form coarse-graining ap- 
proach circumventing the secular approximation): The 
Lamb-shift Hamiltonian for example will only act within 
the subspace of energetically degenerate states. Note 
that when the spectrum of the system Hamiltonian is 
non-degenerate, Eq. (3) may be simplified into a rate 
equation system (which is independent on the Lamb- 
shift) for the diagonals of ps in the energy eigenbasis. 
In the dampening coefficients, the functions 



7a/3 (w) = / C afj (T)e + ^ T dT, 



CTa/3 M 



C^(r)sgn(r) e+^ T dr 



(5) 



are even (jap) and odd (cr a p) Fourier transforms of the 
bath correlation functions 



C aP (r) = Tr B {e +iTH »B a e- irH »B p pB} 



(6) 



The bath correlation functions have many interesting an- 
alytic properties [3]. For example, when the bath is held 
at a thermal equilibrium state (canonical ensemble) 
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one can easily verify [3] the Kubo-Martin-Schwinger [17- 
19] (KMS) condition C a p(r) = C 0a (-T - ifi). Since the 
bath correlation functions are analytic in the lower com- 
plex half plane, the Fourier transform of the KMS con- 
dition reads 



lap{—u) = e /3 "7/3a(+w) , 



(8) 



and can be used to prove [9] that the equilibrated Gibbs 
state 
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Tr { e -/3ffs} 
is a stationary state of Eq. (3). 

II. CONSERVED QUANTITIES 

Now assume that there exists a conserved quantity 
N = Ns + Nb, where N$ and Nb act only on system and 
bath, respectively. That is, we assume that [H$,Ns] = 0, 
[H B ,N B ] = 0, and [H SB ,N S + N B ] = 0, such that non- 
trivial evolution arises via [i?sB,^Ys] 7^0^ [H$b,Nb\- 
The conservation laws imply the identity 
H SB = e~ KiNs+N ^H SB e + ^ Ns+NB K Acting on this 



expression with e 
identity 
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yields with Eq. (2) the 
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such that 
i-kJVs[ ] e 



e 

tor is mapped 
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effectively 



a transformation of the form 
on a system coupling opera- 
into a transformation of the form 
e '"'"[. . .je on the corresponding bath coupling 

operator. We will show in the following that this 
identity leads to additional properties of the dampening 
coefficients in Eq. (3), when the bath density matrix is 
assumed to be in the grand-canonical equilibrium state 
with chemical potential p 



-/3(H B -fiN-B) 
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Note that - depending on the spectrum of H B - normal- 
izability of ps may impose constraints on the chemical 
potential, compare Sec. IV B, Sec. IV C, and Sec. IV D. 

Evidently, as [ifs,JVjy] = 0, we may and will in the 
following choose \a) to be the common eigenbasis of the 
two operators with H$ \a) = E a \a) and N$ \a) = N a \a). 

When we multiply the Lamb-shift coefficients &ij by a 
factor of the form e + P^ Ni ~ Nj \ we may use Eq. (4) to re- 
place the eigenvalues by operators, such that the system 
operators in Eq. (4) are rotated. Then, the identity (10) 
with Eq. (5) can be used to transfer the pseudo-rotation 
to the bath correlation functions. Finally, we may use 
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the invariance of the trace over the bath degrees of free- 
dom under cyclic permutations and [-/Vb,Pb] = to see 
that 



>%3 , 
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i.e., the Lamb shift Hamiltonian only acts on states with 
both degenerate energy and particle number. An analog 
calculation for the dissipativc coefficients ^ a b,cd reveals 
the identity 
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When we consider additional thermal Boltzmann factors, 
one needs to change the integration path in the Fourier 
transform (5) - using that the bath correlation functions 
are analytic - to show that the balance relation 
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holds. Such relations are termed fluctuation theo- 
rems [20]. Specifically, relations (12), (13), and (14) 
generalize the KMS condition (/x = 0) for the quantum 
master equation in Eq. (8) to systems with a conserved 
quantity. 



III. STATIONARY STATE 

A. Single Reservoir 

The matrix elements of the grand-canonical Gibbs 
state (11) read 

Pij — 7? — ^ i (15) 

where Z = Tr je~' 3( -^ fs- ' tArs - ) } denotes the normalization. 
For such a diagonal density matrix, the time-evolution 
of off-diagonal matrix elements may in principle still be 
influenced by the diagonals, as Eq. (3) reduces to 

Pij = -i<Jij (Pjj - pii) + 22 7ia,jaPaa 

a 

a 

To show stationarity of the Gibbs state (15), it is conve- 
nient to distinguish different cases: 

a. Trivially, when i ^ j and also Ei ^ Ej in Eq. (16), 
we have pij = 0, since all coefficients simply vanish, 
cf. Eq. (4). 

b. When i = j and evidently also Ei = Ej, Eq. (16) 
reduces to the rate equation system 



c. For degenerate energy levels we have the additional 
possibility that i ^ j but Ei = Ej . Cancellation of 
the Lamb-shift terms in (16) results from Eq. (12). 
Showing that the remaining dissipative terms also 
vanish amounts to 
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which is directly evident from relations (13) 
and (14). 

To summarize, we have shown that the state (15) is a sta- 
tionary state of the quantum master equation (3), when 
the reservoir density matrix is of the form (11). Gener- 
ally of course, the existence of further stationary states is 
possible, but for an ergodic [3] evolution the BMS approx- 
imation scheme for a single reservoir leads to equilibration 
of both temperature and chemical potential. This equili- 
bration has been noted earlier for specific examples [21- 
23] and has been used quite generally for systems of rate 
equations [24]. Here however, we have a rigorous proof 
for the quantum master equation. 



B. Multiple Reservoirs 

When the system of interest is is not only coupled to 
a single, but multiple (K) reservoirs 
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where varying coupling strengths are absorbed in the 
Ba operators and the independent reservoirs are char- 
acterized by different inverse temperatures (3^ and dif- 
ferent chemical potentials fj,^ k ' 
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much less is known about the resulting stationary 
state [25]. A decomposition of the interaction Hamil- 
tonian in the form of Eq. (19) with identical system 
coupling operators for each bath is always possible, as 
we have chosen the A a operators to form a complete 
basis set for hermitian operators in the system Hilbcrt 
space. We assume that some interaction Hamiltonians 
may obey a conserved quantity = N$ 



Ng , where 



(fc) rr(k) 
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Jia.iaPaa — ^ Jai,aiPii = , (17) dently, the form of Eq. (3) remains invariant with 
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which vanishes due to the detailed balance rela- 
tion (14), evaluated for i = j. 
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where 7^ cd and describe the dissipation and Lamb- 
shift, respectively, due to the fc-th bath only. Accord- 
ingly, each bath yields separate detailed balance condi- 
tions of the form of Eqns. (12), (13), and (14). In general, 
this will lead to a non-equilibrium stationary state. 



C. Rate Equations for Ladder Spectra 

Quite general statements on the resulting non- 
equilibrium stationary state may be obtained for mas- 
ter equations that assume tri-diagonal form in an 
A-dimensional energy and number cigenbasis (with 
Hs\m) = E m | to) and N s \m) = N m |m)), where p m = 
(to | ps | to) and we assume ordering with respect to the 
(quasi-)particlc numbers in the system (N m+ i — N m = 1) 

Pm ^m.rn—lPm—1 ^m.m+lPrn+l 

[7m-l,m "i" Tm-(-l,m] Pm • 

(22) 

In this basis, the populations evolve independently from 
the coherences (which we assume to decay), and tran- 
sitions between populations arc mediated by single par- 
ticle tunneling processes. Note however, that even an 
effective rate equation system of the form (22) may 
keep genuine quantum properties as the eigenstates |to) 
themselves may e.g. be entangled between different sub- 
parts of the system. At the boundaries toi and of 
the spectrum (where N mi — or N m2 = N) the un- 
physical tunneling rates vanish, for example we have 
p mi = +lm u m 1 +\Pm l +i-lm l +i,m l Pm l - Computing the 
stationary state and using the results from the boundary, 
this implies that for all to, it has to satisfy 



Pm + l _ 7m+l,ra 
Pm ^m^m+l 



(23) 



Together with the trace condition P~ m = this com ~ 
plctcly defines the stationary state. In addition, we as- 
sume that (compare Sec. IV for examples) the dampening 
coefficients associated with the k th bath have the decom- 
position 
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where 



E m+1 - E m , and g m , G (fc) (w) and 

F±\w) contain the details of the system, the cou- 
pling, and the thermal properties of the bath, re- 
spectively Using this factorization in the local bal- 
ance condition (14) leads with N m+ i — N m = 1 to 

/m,m+l c 

matically fulfilled for fermionic 
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which is auto- 
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f (fc) H (25) 



or bosonic 



(w) (26) 



baths. It is obvious from Eq. (23) that for coupling to 
a single bath, neither g m nor the energy dependence of 
the tunneling rate (w) do affect the stationary state 
- only the transient relaxation dynamics will be changed. 
However, when the system is coupled to multiple baths 
obeying Eqns. (24), it is easy to show that the stationary 
state - as characterized by Eq. (23) 
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1 ± F±(w m+ i, m ) 



(27) 



is the same as one for a hypothetical single structured 
bath with the weighted average occupation function 



(28) 



For a single (K = 1) bath of either fermionic or bosonic 
nature the ratio in (27) reduces to the conventional Boltz- 
mann factor p m +i/pm = ("m+i.m-M ^ sucn that 

we will term it generalized Boltzmann factor further- 
on. Evidently, the energy dependence of the tunnel- 
ing rates enters the average occupation function (28) 
and may therefore be used to tune the resulting non- 
equilibrium stationary state. For example, whereas the 
canonical Gibbs state (9) will for finite temperature al- 
ways favor the ground state and the grand-canonical 
Gibbs state (11) may favor states with a certain particle 
number, it is here possible e.g. to select multiple states 
of interest. Note however, that with using only bosonic 
baths (with n(w) > 0) it is not possible to achieve e.g. 
pm+i > Pm, which is in stark contrast to fermionic baths, 
where this only requires F_ (w m +i. m ) > 1/2. 

To illustrate this idea let us for simplicity parameterize 
the tunneling rates phenomenologically (but see e.g. [26, 
27] for a microscopic justification) by a Lorentzian shape 



r<*>(w) = 



(w - 0o k ) 2 + S 2 k 



(29) 



with maximum rate at frequency uj^ and width 5k- 
When the average occupation function at the transi- 
tion frequency ui m+ i^ m between states \m) and |to+ 1) 
exceeds 1/2, the population of the states will increase 
Pm+i > Pm, whereas the opposite is true for i{io m +i.m) < 
1/2, see Fig. 1. Depending on the number and the 
thermal properties used baths, the resulting hypotheti- 
cal average distribution (28) may assume quite arbitrary 
shapes, such that more sophisticated statistical station- 
ary mixtures than in Fig. 1 are conceivable. 
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FIG. 1: (Color Online) Effective bath occupation (solid black) 
in Eq. (28) generated by two fermionic baths (dotted red and 
dashed green curves, chemical potentials // fc ' and temper- 
atures ^ fc 'J are given by the center and width of the 

shaded boxes, respectively) with Lorentzian tunneling rates 
of form (29). The detuning of the tunneling rates is ad- 
justed such that the average hypothetical occupation func- 
tion f(aj) is not monotonically decaying. Then, the thresh- 
old 1/2 is passed several times (vertical solid blue lines) 
within the transition frequency window probed by the system. 
This will lead to an increasing population p m +i > pm when 
f(w m +i, m ) > 1/2 and to a decreasing population p m +i < pm 
when f(a; m +i jm ) < 1/2 (compare the inset). Such multi- 
modal distributions are clearly non-thermal, i.e., they cannot 
be generated by a single thermal bath by adjusting tempera- 
ture and chemical potential. Parameters have been chosen as 
Ae £ (1) = Ae £ (2) = 2, /x« = Ae 0) p (2) = 9Ae , I\ = T 2 , 
5i = 82 = O.lAeo, and o>i = Aeo with u)2 = 9Aeo- The shown 
equidistant transition frequencies (vertical solid and dashed 
blue lines) correspond to eigenvalues that scale quadratically 
with m and may be realized by the example in Sec. IV A by 
using N = 10, e = Ae , U = Ae , and T — 0. 



fix the parameters /3 and p. Note however, that then in 
general J3 < is possible, which corresponds to popula- 
tion inversion (i.e., favoring for p, = the most excited 
state) Evidently, this forbids the interpretation of the 
parameter (3 as inverse temperature. 

Alternatively, for N > 2, Eqns. (27) may still be lin- 
early dependent, as is the case when the system Hamil- 
tonian provides only a single transition frequency (e.g., a 
harmonic oscillator) . In an approximate sense, linear de- 
pendence may also be generated when the average occu- 
pation function F±(to) is essentially flat in the transition 
frequency window probed by the system. Assuming only 
a single transition frequency il, such that the tunneling 
rates are described by a single number G^ k \fl) = Tk, 
and vanishing chemical potentials pS k ^ = = p, in the 
low-energy limit fi^Vl <C 1, Eqn. (27) defines average 
temperatures for bosons and fcrmions 
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E T i Tik) . 
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r T( fc ) 



(31) 



with r = y\ Ti, that correspond to the weighted arith- 
metic (bosons) and weighted harmonic (fermions) mean. 
The bosonic average temperature does well resemble the 
Richmann mixing formula and has been found previ- 
ously [10]. However, it should be noted that this is 
different from the so-called classical (high-energy) limit 
@W fl > 1, for which one obtains an arithmetic mean of 
the Boltzmann factors from Eq. (27) 



(32) 



for - as expected - both fermionic and bosonic baths. 



D. Trivial dynamic equilibria 

It has been observed for special systems [10] that the 
equilibrium state had a thermal form with some average 
temperature even though the system of interest was cou- 
pled to more than just a single thermal bath with differ- 
ent temperatures. For the ladder-like systems discussed 
in the previous section, it is obvious that to describe the 
stationary state with just two parameters j3 and /2 re- 
quires Eqns. (27) to be compatible with 

Pm 1 ± F± (uj m+l,m ) 

where the effective average occupation is defined in 
Eq. (28). This can be achieved in several possible ways: 
Firstly, for a three-dimensional system Hilbert space 
(N = 2 in the previous section), we will only have two 
transition frequencies and accordingly only two equations 
of the above form. These two equations would uniquely 



IV. EXAMPLES 

For a single reservoir the equilibration of both tem- 
perature and chemical potential has been observed for 
interacting double dots [21] in the BMS approximation. 
Therefore, we only give examples to illustrate the results 
in Sec. IIIB, Sec. IIIC, and Sec. HID. Naturally, in case 
of only a single coupling bath, the case of Sec. Ill A is 
also reproduced here. 



A. Homogeneous Electronic Nanostructure 

Consider a nanostructure with N homogenous elec- 
tronic sites (we neglect the spin) 

N U 

H S = e E 4di + -E 4*$dj + T E 4<h ( 33 ) 

i=l i^j i^j 
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with single-particle energy e, Coulomb interaction U, 
and hopping term T that are all assumed com- 
pletely isotropic. The permutational symmetry sug- 
gests to reduce the problem to the symmetrized subspace 
with the basis \N,m+l) = —, 1 \N, m) with 

y (N—m)(m+l) 

< m < N denoting the number of electrons in the sys- 
tem, where D = ^2 i=1 di and \N,0) = |0, ...,0) repre- 
sents the N particle Fock space vacuum. Clearly, both 
Hs and N$ = J2i d-ldi are diagonal in this basis. For 
N = 1 we recover the single resonant level [28], but for 
for N > 1, the spectrum of the system Hamiltonian be- 
comes nontrivial. The eigenvalues of the symmetric sub- 
space read E m = me + m(m — l)U/2 + m(N — m)T, 
such that the spectrum may become equidistant when 
U = 2T with Lu m +i,m -» e + (N - 1)T - which still ad- 
mits a strongly interacting model. 

At first we assume that the nanostructure is only cou- 
pled to a single lead Hb = Ylk £ fc c ! c fc a ^ temperature 
P and chemical potential fi via the tunneling Hamilto- 
nian [33] H SB = D®Y,k 9kc\ + £> f <8> 9 k °k, where g k 
represents a frequency-dependent coupling constant, and 
Cfc are fcrmionic annihilation operators acting on the lead 
Hilbcrt space. The conserved quantity is composed from 
Ns and Nb = ^2 k c k c k . We may write the interaction 
Hamiltonian also as H$b = A\ ® B\ + A 2 <S) B 2 , with 
the hermitian and trace-orthogonal (not-normalized) 
system coupling operators A\ = (D^ + D) and 
A 2 — i (£)t — D) . and the associated hermitian 

bath coupling operators B\ = (^g k c\ + g k Ckj /2 and 

B2 = i J2k (fffc ! — Sk c kj /2- For a bath in thermal 
equilibrium with inverse temperature /3 and chemi- 
cal potential p, such that its density matrix is given 
by Eq. (11), we obtain for Fourier transforms (5) of 
the bath correlation functions 711 (w) = 722 (w) = 
T{+w) [1 - f(+w)] /4 + r(-w)f(-w)/4 and 7ia (w) = 
7| 1 (u;) = iT(+u) [1 - f(+w)] /4 - iT(-w)f(-u})/4, where 
T(u>) = 2nJ2k \9k\ 2 S(u!~ek) is the tunneling rate and the 
Fermi function f(w) encodes the bath properties /3 and 
p, compare Eq. (25). Obviously, the Fourier-transform 
matrix of bath correlation functions has non-negative 
eigenvalues - a consequence of Bochners theorem [3, 29]. 
These lead to the non- vanishing dampening coefficients 

7m,m+i = (N - m)(m + l)r(w m+ i, m ) [1 - f(w m +i, m )] , 
= (N - m)(m + l)r(w m+ i )fn )f(a; m+ i, m ) , (34) 

where w OT +i im = £, n +i — E m and the factoring condi- 
tion (24) is obviously fulfilled. We obtain a rate equation 
of the form (22) with p m = (N,m\ps \N,m), where the 
Lamb-shift terms are irrelevant, since we have by exploit- 
ing the permutational symmetry mapped our system to 
a nondegenerate one. 

Now we consider tunnel couplings to multiple baths 
with factorizing density matrices as in Eq. (20). The 
form of Eq. (22) remains invariant, and we simply have 

7m, m±i = Z)fc7m,m±D witn different temperatures /3 (fc) 



and chemical potentials /i^ entering the rates as in 
Eq. (34). The general non-equilibrium steady state of 
Eq. (22) fulfills 



Pm+l _ Efc r(fc) ( w m+l,m)f (fe) (w m+ i im ) 



Pn 



Efc T« (W m+ l, m ) [1 - f « (Wm+l,m)] 



(35) 



and is thus identical with the non-equilibrium steady 
state for coupling to a single hypothetical non- 
equilibrium bath with average non-thermal distribution 
of the form (28), compare also Fig. 1. 



B. Coupled oscillators 

We consider a harmonic oscillator Hs = flb^b cou- 
pled to many others Hb = Sfc w *-'^i^fc with positive 
eigenfrequencies ujk > via quasi-particle tunneling 
H SB = b <g> Y, k h k b l + ^ ® Efc K b k- The conserved 
quantity is composed from Ns = b^b and Nb = ^ k b\b k . 
Rewriting the interaction Hamiltonian in terms of her- 
mitian operators, we obtain A\ = {br + b), A 2 = i(b* — b), 
Bi = E k (hk b l+h* k b k )/2, and B 2 = 1 J2 k (h k b{-h* k b k )/2. 
The matrix elements of the Fourier transforms (5) of the 
bath correlation functions equate to 711 (w) = 722 (w) = 
1/4 [+e(+w)[i + n(+w)]r(+w) + 0(-w)n(-w)r(-w)] 
for the diagonals and 7^ (w) = 712 (^) = 
i/4 [+0(+w)[l + n(+w)]r(+w) - 0(-w)n(-w)r(-w)] 
for the off-diagonals, where &(oj) denotes the Heaviside 
step function, T(u>) = 27r^ fe (5(cj — uj k ) the quasi- 
particle tunneling rate, and n(cj) denotes the bosonic 
occupation number as defined in Eq. (26). The condition 
that p < minfe(a;/ s ) grants positivity of all bath occupa- 
tions. In addition, it implies that r(f2)| n< = 0, such 
that we assume also fl > p throughout. Accordingly, 
the Fourier transform matrix of the bath correlation 
functions is positive semidefinite at all frequencies. In 
the Fock space basis (where b' \n) = \Jn + 1 \n + 1)), 
we obtain a rate equation of the form (22), where 
< n < 00. However, since even the original cigenstates 
arc non-degenerate, the dampening coefficients equate 
(with Q > 0) to 



7„,„ +1 = (7i + l)r(fi) [1 + n(fi)] , 
J n +i,n = (n + l)r(fi)n(fi) , 



(36) 



which is also compatible with assumption (24). The re- 
sulting system is infinitely large, one may however, intro- 
duce a cutoff size N cut and solve for the stationary state 
of the rate equations for finite A^ cu t. The ratio p n +i/p n of 
two successive populations yields the desired Boltzmann 
factor, which even happens to be independent of N cut . 

For multiple baths, we obtain equilibration in a ther- 
mal state with the unique generalized Boltzmann factor 



Pn 



X; fc rW(Q) n W(») 
~ ~ £ fc r«(ft) [l + n ( k )(ft)] 



(37) 
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consistent with Eq. (27). This generalized Boltzmann 
factor is the same that one would obtain for contact with 
a single hypothetical non-thermal bath at an average oc- 
cupation compatible with Eq. (28), which in the high- 
temperature and p^ — > limit reduces to the bosonic 
average temperature in Eq. (31). This coincides well with 
the well-known temperature mixing formula and previous 
results [10]. 



C. Spin-Boson Model 

A variant of the spin-boson model coupled to a sin- 
gle bath has already been provided in Rcf. [30], such 
that we here only generalize to K > 2 baths and non- 
vanishing chemical potentials. We consider a large spin 
system H§ = VL/2J Z with ft > 0, coupled to a bath 
of harmonic oscillators Hb = J2k w k b\b k with 0J k > 
via H SB = J + ® Y,k hkb' k + J~ ® J2k hkbk, where J a = 
J2iLi a ? and cr± = (c x i i& v )/2. The conserved quan- 
tity is given by N s + N B = -J z /2 + *£ k b\b k . We im- 
pose the same conditions on the chemical potential(s) 
as before: p < mhi k (u} k ) and fx < 0. Choosing the 
system coupling operators as A\ = J x and A2 = J y . 
the Fourier transforms (5) of the bath correlation func- 
tion are identical to the previous section. In order 
to calculate the dampening coefficients, permutational 
symmetry suggests to use the angular momentum ba- 
sis \N/2,m) with -N/2 < m < +N/2. Using that 
\ j,m) = + 1) — m(m ± 1) \j,m± 1), we obtain 

a rate equation of the form (22) with the coefficients 



^m+l^m 



N ( N 



2 V 2 
N ( N 
~2 { ~2 



- m(m + 1) 
+ 1 I - m(m + 1) 



r(n)n(fi) . (38) 



Solving this rate equation with multiple baths for its 
steady state yields a thermal bath with the same unique 
generalized Boltzmann factor as Eq. (37), as predicted 
by Eq. (27). 



D. Mixed Spin Model 

We consider a spin- 1/2 system H§ = VL/2a z that is 

firstly coupled to a bosonic bath H-£' = J2k u kb k b k 

via the dissipative coupling fflg = <r + J^k ^kb k + 

o~ g) J2 k h* k bk, and secondly to a fermionic bath H B = 

J2k e kc{ c k via the coupling H<£ — a + (g> Y^k 9kc\ + 
a~ ® 12k9k c k- Note that in contrast to the previ- 
ous examples we do now consider two different baths 
from the beginning. The interaction Hamiltonians ex- 
plicitly obeys the conserved quantity constructed from 

Ns = -<r*/2, = Ekblbk, and = EA<*- 

(2) 

Note that i?g B docs not conserve the number of fermions, 



but such a model may represent scattering processes 
with a further fermionic bath that omitted from the 
description. As before, we require that ^S 1 ' < f2 and 
< mm k (ui k ). Choosing the system coupling op- 



/' 



(i) 



erators as A\ = a x and A2 = u v , we obtain B\ 
l/ZEJhkbl + htbk), B 



(i) 



^Ek[hkbl-hlbk 



and b[ 2) = l/2J2 k [9kcl+g* k c k ) with B 



? (2) 

'2 



i/2 J^k [dkcj, — g k c k ) ■ The Fourier transform of the bath 
correlation function for the first bath corresponds to 
Sec. IV C, whereas the Fourier transform matrix for the 
second bath is identical to that of Sec. IV A. Accord- 
ingly, we obtain in the cr z -eigenbasis a z \a) = (— l) a \a) 
with a € {0, 1} and p a = (a\ ps \a) the master equation 



Po = - 




'1,0 


(2)' 

1-7$ 


Po + 




h7@' 


Pi = - 


f- 


V 1} - 

'1,0 


r 71,0 


Po - 


7g - 





Pi 



Pi 



(39) 



with the dampening coefficients 7^ = ri(O) [1 + n(f2)], 

7$ = ri(n)n(n), 7$ - r 3 (n)[i-f(n)], and 7$ = 

r2(Q)f(f2), where f(w) and n(w) have been defined 
in Eqns. (25) and (26), respectively, and T\{us) = 

^Y,k\ h k\ S(oj- io k ) with r 2 (w) = 2irJ2 k \9k\ 5(oj - e k ) 
represent the coupling strengths to the two baths, respec- 
tively. The stationary state of Eq. (39) is characterized 
by the generalized Boltzmann factor 



Pi 



ri(n)n(n) + r 3 (fi)f(n) 



Po ri(fi)[i + n(n)] + r 2 (n)[i-f(n)] 

which is consistent with Eq. (27). 



(40) 



V. CONCLUSIONS 

Under the Born, Markov, and secular approximations, 
quantum systems coupled to a single bath described by 
inverse temperature /3 and chemical potential p relax - 
when the total Hamiltonian conserves a (quasi-)particle 
number - into a stationary equilibrium ensemble that 
is described by the same inverse temperature and the 
same chemical potential. As long as only a single bath is 
involved, this also holds when the spectrum of the system 
Hamiltonian is not equidistant. 

For coupling to multiple thermal baths and tridiago- 
nal rate equations, a hypothetical non-thermal average 
bath is effectively felt by the system, which will in gen- 
eral lead to a non-thermal stationary state. However, 
when there exists only a single transition frequency as 
e.g. in two-level systems, the resulting stationary state 
may be well characterized by a single Boltzmann factor 
with two parameters (3 and p. The possibility of creating 
level inversion however demonstrates that then /? does 
not always define an inverse temperature. 

There are several interesting consequences: Firstly, 
by using equilibration under side constraints one should 
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be able to prepare not only the ground state of a sys- 
tem Hamiltonian by dissipative means but also an en- 
ergetically sufficiently isolated energy eigenstate with a 
desired particle number when temperature and chemi- 
cal potential of a single grand-canonical bath are tuned 
accordingly. Secondly, in order to generate interesting 
non-equilibrium stationary states via coupling to multi- 
ple baths, fermionic baths appear more favorable. Be- 
yond this, it is necessary (though not sufficient) to con- 
sider system Hamiltonians with multiple allowed transi- 
tion frequencies. Energy-dependent tunneling rates can 
then significantly enhance the non-thermal signatures of 
the resulting stationary state. Finally, when the focus 
is on particle or thermal transport, the quantities of in- 



terest are often the stationary current through the sys- 
tem and its fluctuations (noise), see e.g. Refs. [31, 32]. 
The calculation of these quantities heavily depends on 
the knowledge of the stationary state and may therefore 
strongly benefit from its analytic knowledge. 
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